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ABSTRACT 

We study the orbital evolution of wide binary stars in the solar neighborhood 
due to gravitational perturbations from passing stars. We include the effects 
of the Galactic tidal field and continue to follow the stars after they become 
unbound. For a wide variety of initial semi-major axes and formation times, 
we find that the number density (stars per unit logarithmic interval in projected 
separation) exhibits a minimum at a few times the Jacobi radius rj, which equals 
1.7 pc for a binary of solar-mass stars. The density peak interior to this minimum 
arises from the primordial distribution of bound binaries, and the exterior density, 
which peaks at ~ 100-300 pc separation, arises from formerly bound binaries that 
are slowly drifting apart. The exterior peak gives rise to a significant long-range 
correlation in the positions and velocities of disk stars that should be detectable 
in large astrometric surveys such as GAIA that can measure accurate three- 
dimensional distances and velocities. 

Subject headings: binaries: general — Galaxy: kinematics and dynamics — solar 
neighborhood — stars: kinematics 



1. Introduction 



Wide binary stars are disrupted by gravitational encounters with pa ssing s tars, m olec- 
ular clouds, and other perturbers. This process was first investigated by lOpiki (119321 ). who 
estimated the e-folding time for disruption of binaries composed of solar-mass stars, with 
apocenter distances of 1 pc , to be 10 Gyr or less. O t her ea r ly estimate s of th e disrupt i on rat e 
are due to lAmbartsumianl (119371 ) , IChandrasekharl (119441 ) , lYabushital (119661 ) , iHeggid (119751 ) , 
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Kind fll977h . lHeggid (119771 ) . iRetterer fc Kind (119821 ) . and lBahcall. Hut, fc Tremaind (119851 ) . 
In the last of these a binary of age t and component masses Mi and M 2 was estimated to 
have a 50% survival probability at semi-major axis 



h/2 



(to) = 0.002 



(Mi + M 2 )a 



3.1 x 10 4 AU 



M 1 + M 2 a 0.03M| pc" 3 10 Gyr 



Gp 2 t "' 2M 50 km s ' p 2 

Here 3a 2 is the mean-square relative velocity between the center of mass of the binary and the 
perturbing stars, and p 2 = J n(M p )M 2 dM p is the second moment over mass of the number 
density of stars in the solar neighborhood (cf. eq. I4"7l) . This estimate is substantially shorter 
than Opik's, mostly because it includes the cumulative effects of distant, weak encounters 
(which Opik recognized to be important but did not compute). 

A closely related problem is to estimate the distribution of semi-major axes of wide 
binaries. For relatively small semi-major axes a < ^1/2(^0) the distribution is presumably 
primordial, and thus reflects the (poorly understood) formation process of wide binaries. At 
larger semi-major axes a > ai/ 2 (t ) the distribution ought to be primarily determined by the 
disruption process. The Fokker-Planck equation that descri bes the evo lu tion of the semi- 
major axis distribution (eg . 1351) was derived and solved by iKind (119771 ). IRetterer &: King 
(119821 ) and IWeinberg et ajj (Il987l ). who showed that dn oc da/a 2 for a > ai/ 2 (to)- 



tr 



(1) 



Observationally, the distribution of wide binary semi-major axes is determined by mea- 



suring the projected separations of common-proper-motion binaries 



e.g. ; 



Chaname fc Gould 



2004 . Poveda et al. 2007 . Lepine fc Bongiorno 2007 . Sesar et al. 2008 ; see also Chanami 



2007 



and references therein). For a < 3 x 10 3 AU the distribution of separations or semi-major 
axes3 of disk binaries is approximated well by Opik's (1924) law, 

da 



dn oc d log a 



(2) 



At larger semi- major axes, the nu mber of binaries falls more s teeply, roughly as dn oc da /a 



1.6 



for 3 x 10 3 AU < a < 10 5 AU ( ILepine fc Bongiorno! 120071 ) . The further steepening to 
dn oc da I a 2 that is expected for a > 
cult to detect. There have been a nu 



'1/2 



10 Gyr) 



3 x 10 4 AU is much more dim- 



w^v,„~, ^ ^ 

Wasserman & Weinberd 


1987 


Latham et al 


|l991 


Yoo. Chaname & Gouldl 


2004 


Quinn et al. 


2009) 



are controversial ( 


Bahca 


& Soneira 


1981; 


Wasserman & Weinberg 


1991 


Palas 


2000 



tribution are likely to improve dramatically in the next few years because of large, accurate 



For a population of binaries at a given semi-major axis a, with other orbital elements assigned as 
described at the start of $31 the median projected separation is 0.978a. Thus we may assume that the 
distributions of semi-major axes and separations are nearly the same. 
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proper-motion surveys. In particular, the GAIA spacecraft will determine both proper mo- 
tions and trigonometric parallaxes for millions of nearby stars with unprecedented accuracy, 
allowing a far better determination of the binary population at large separations than the 
ground-based proper motions and photometric parallaxes that have been used in all studies 
so far. 



Large, well-characterized samples of wide binaries have many applications (jChaname 



20071 ). In particular, the distribution of wide binary semi-major axes can be used to con- 
strain the properties of molecular clouds and other mas sive structures in the disk, and 



possible compact objects (M ACHOs) in the dark halo (IBahcall. Hut. &: Tremaind 11985 



Yoo. Chaname fc Gouldll2004l ). If the distribution of binaries can be measured at separations 
as large as a few parsecs we expect to see "tidal tails" of the kind that have been detected 
aroun d globular clusters ( Odenkirchen et aLlboOl ; Belokurov et al.|[2006 ; Grillmair &: Dionatos 
20061 ); the evolution of these structures offers a prototype for the evolution of the phase- 
space structures in the solar neighborhood caused by the disruption of stellar clusters 
JPehnen fc Binnevlll998h . 



Almost all theoretical studies of the expected distribution of wide binaries have made 
two related approximations that compromise their validity at the largest semi-major axes: 



• The stars are assumed to disappear instantaneously as soon as their orbits become 
unbound. This is unrealistic because the disruption rate is dominated by weak, distant 
encounters, so most escaping stars have very small relative velocity and only drift 
slowly apart. 

• The Galactic tidal field is ignored. The tidal field becomes stronger than the gravita- 
tional attraction between the stars in the binary when the separation is roughly the 
Jacobi or tidal radius, which equals 1.7 pc = 3.5 x 10 5 AU for solar-mass stars in the 
solar neighborhood (eq. [43]). Thus the tidal field is already significant at the separa- 
tions (~ 10 5 AU) probed by current measurements of the wide binary distribution, and 
dominates the dynamics at larger separations. 



Including these two effects is necessary if we are to understand the expected distribution 
of binary stars — bound and unbound — at semi-major axes of 10 4 AU and larger. To achieve 
this understanding is the primary goal of this paper. We restrict ourselves to the evolution 
of disk binaries under the influence of passing stars, although it is straightforward to extend 
our methods to include either halo binaries or other perturbers such as molecular clouds or 
massive black holes. 
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The structure of this paper is as follows. In we describe the basic equations of motion 
for binary stars in the Galactic tidal field and how we calculate the perturbations from other 
stars that drive the orbital evolution. We also review the standard analytic treatment of the 
evolution of bound binaries using a diffusion equation. Then in §3] we describe the results 
from our simulations. Finally, §H] contains a discussion and conclusions, and Appendix [A] 
derives an analytic model that approximately describes the diffusion of unbound binary stars. 



2. Basic equations in the numerical simulation 

In this section, we describe the details of our numerical simulation. First, we give the 
equations of motion of the binary star in Hill's approximation. Then, we describe how we 
include the effects of kicks from other stars. Finally, we describe the diffusion approximation, 
which should be valid for binaries with small semi-major axes. 



2.1. Evolution without kicks 



We use Hill's approximation (e.g.. iHeggid l200ll ; iBinney &: Tremaind 120081 ) to describe 
the motion of the binary star in the Galaxy. Hill's approximation is valid because the mass 
of the binary is much less than the mass of the Galaxy (by a factor ~ 10 11 ). Let the masses 
of the two stars in the binary be Mi and M 2 . We assume that the potential of the Galaxy 
is symmetric about the plane Z = 0, where (X, Y, Z) or (R, 0, z) is an inertial Cartesian 
or cylindrical coordinate system with origin at the center of the Galaxy. We introduce a 
second coordinate system (x, y, z) with origin in the Z = plane at distance R g from the 
Galactic center. The x-y and X-Y planes coincide but the origin of the (x, y, z) coordinate 
system co-rotates with the Galaxy. The x-axis points radially outward, the y-axis points in 
the direction of Galactic rotation, and the z-axis is perpendicular to the Galactic plane. The 
x, y and z axes form a right-hand coordinate system. With these conventions, the positive 
z-axis points toward the South Galactic Pole. The angular speed of the Galaxy at radius 
R g , which equals the angular speed of the (x, y, z) frame, is Q g = Q g e z in the z or "vertical" 
direction. The angular speed fl g is related to the potential of the Galaxy $o(R, z) by 



ft 2 



R„ 



(3) 



where <&' (R g ,0) = d&o/dR\m gt Q). As we want to study binary stars in the solar neigh- 
borhood, we can just choose R g to be the distance of the Sun from the Galactic center, 
R g = 8 kpc. 
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In the co-rotating frame with origin at R g , the position of star i, % — 1, 2, is labeled by 
r-j = (xi,yi, Zi). Then in the co-rotating frame with origin at the center of the Galaxy, the 
star's position is H g + r^, where R 9 = (R g ,0, 0) and the equation of motion for either star 
in the binary system is 



d 2 (ri + R 



2f2 9 x 



tZ(r< + R s 



o 9 x [n fl x (n + r s )] , 



(4) 



dt 2 ' °~ ~" 9 ' dt 

where Vj is the gradient with respect to r-j. The potential $ includes the contribution from 
the Galaxy $o a s well as the potential of the binary stars For <£>o, we use the distant- 
tide approximation, which means §o(R g + x, y, z) at the position of a star is expanded with 
respect to $ g (R g , 0, z). Then we have 

(V$ ) a = $o, a + ®o,ap i3 + 0(r 2 ), a=(x,y). (5) 

P=x,y,z 

Here $ ,a = (d®o/d a )(R g ,o) an d §o, a /3 = (<9 2 $o/<9 a dP)(R g ,o)- Equation is correct for any 
value of r^; in particular, when = 0, it is correct for the center R 9 , and we may subtract 
the equation for the center from (OH) to obtain 

d 2 x d$ b 
~d¥ = 



dt 2 

cPz_ 

dt 2 



dx 
dy 



(3=x,y,z 



„ dr 

2Cla X — 

9 dt 



- J2 ®°>ypP 

P=x,y,Z 

d$ 



„ dr 

x — 

9 dt 



We have 



dz dz 



/<9 2 $ 



2S7„ x — 

9 dt 



- x x r)] a 

" x x r )] 4 
fi 9 x (O s x r)] 2 . 



J y 



V Si? 2 J (Rg,o)' ^°' yv \R dR J (Rg,0) 
The potential for star % — 1, 2 is just the potential from the other star in the binary. Then 
we have 

GM 2 



1 <9$ r 



0,xy 



0. 



(6) 



(7) 



61 



a/ (xi - x 2 ) 2 + (yi - y 2 ) 2 + (zi - z 2 ) 2 
where (xi,yi, Zi) denotes the position of star i, and the formula for $ b2 is obtained by inter- 
changing 1 and 2 in all subscripts. Then the equations of motion for either star are 



x\ =2tt g y~i + 



^-$0(^,0) 



dx,. 



dzi 



Q 2 



Rn 



dy % 



<9$ 

dz,. 



bi 



(9) 
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As the angular speed Q g is related to the potential $o via equation we have 



dtt 



%(R g ,o) = n 2 g + 2R g n g — 



As usual, the Oort constant A(R) is defined as 

We label A g = A(R g ). Then equation (jHJ) can be simplified to 

dyi 



(10) 



(11) 



£i + 



<9$0 



dZi 



From equation (151) , we have the following relations 



Mr 



-M 2 — ^, My 



ox\ 0x2 oyi 01/2 oz\ 

The center of mass of the binary system r cm is defined to be 

M 1 r l + M 2 r 2 



= —M- 



<9$ 



62 



(12) 



(13) 



1 cm 



M 1 + M 2 ' ^ 14) 

The relative coordinates of the two stars are r = ri — r 2 . By adding equations (j!2p multiplied 
by appropriate coefficients for star i = 1 and i = 2, we get the equation for the motion of 
the center of mass^l 

jjcm 2figA C m j 

<9$ 



^cm "I - n 



=0 . 



(15) 



Due to the symmetry of the Galactic potential, (d§o/dz) z= o = 0, so for stars not very far 
from the mid-plane of the Galaxy, we approximately have d$o/dz = (d 2 &o/dz 2 ) z=Q z and 
we define 

(16) 



dz 2 



(R g ,0) 



2 We assume the separation of the two stars along the z direction is much smaller than the thickness of 
the Galaxy. 
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where v g is the frequency for small oscillations in z. The general solution to the above 
equations of motion for the center of mass is just epicycle motion, 

x cm {t) =x 9tCm + Xcos(K g t + a) , 

2Q 

y C m{t) =y g ,cm(t) - Y sin(K g t + a), y g , cm {t) = y gfl - 2A g x 9tCm t, Y = — 9 -X, 

Kg 

z cm (t) =Z cos(u g t + a z ). (17) 

Here Xg^cmi X j yg :0 , Z, a, a z are arbitrary constants, and k 9 is the epicycle frequency defined 
by 

«g = 4fl,(n tf - 4,) . (is) 

The variables 

%g,cm(t) and Ug t cm 

(t) give the position of the guiding center — the center of the 
epicyclic motion — for the center of mass. Subtract equations (1121) with i = 2 from i = 1 and 
we get the equations for the relative motion of the two stars 



x — 2Q g y — AVLgAgX = — 
y + 2Q g x = — 



G(M 1 + M 2 )x 
[x 2 + y 2 + z 2 ) 3 / 2 

GjMi + M 2 )y 
[x 2 + y 2 + z 2 ) 3 / 2 



z | U 2 Z _ G(M 1 + M 2 )z 

+ B ( X 2 +y 2 + Z 2 ) 3 / 2 ■ 1 ' 

In this set of equations, the terms involving y or x arise from the Coriolis force, as we are 
working in a rotating frame; the terms involving A g or v g represent the effect of the Galactic 
tide, and the terms on the right side represent the gravitational force between the members 
of the binary system. 

Equations (fl9l) show that the relative motion is the same as that of a test particle around 
an object with the mass M\ + M 2 in the Galactic tidal field. A special solution to the above 
equations is the stationary solution (x = x = y = y = z = z = 0) 



y = z = 0, x = ±rj, where r,j = 



G(M 1+ M 2) r 



4{l,A g 



is the Jacobi or tidal radius of the binary system. The stationary points are actually the 
Lagrange points in the three-body system composed of binary star and the Galaxy. As will 
be seen from our simulations below, the Jacobi radius sets the characteristic scale for the 
distribution of binary stars at large radii. 
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Equations (fT9l) admit one integral of motion, the Jacobi constant 

p J/.J, .2, .J An A 2 ± 2 2, G(M l + M 2 



= ^(i 2 + 2/ 2 + ^) + $cff(x,i/,z), (21) 

where <& eS (x, y, z) = v 2 z 2 /2-2Vt g A g x 2 '—G(M 1 +M 2 ) / ' \Jx 2 + y 2 + z 2 is the effective potential. 
The Jacobi constant for the stationary solution (120]) is called the critical Jacobi constant E c , 
and is given by 

E c = $ cfr (±0, 0, 0) = -2Q g A g r 2 - G{Ml + M ^ = _ ^-(n gAg ) 1 / 3 [G(M 1 + M 2 )} 2 '\ (22) 

As x 2 ,y 2 ,z 2 > 0, the motion is constrained to the region in which $ e fj (x, y, z) < Ej, and 
the boundary of this region is the zero- velocity surface for a given Jacobi constant, defined 
implicitly by 

$ ea (x,y,z)=E J . (23) 

We choose the time unit to be l/f2 s and the length unit to be rj. Then we can define 
the following dimensionless variables (see eq. 0H for numerical values of the scaling factors) 



r r 

1 -/ 1 ~n 



r" = — — . (24) 



AA g 


X 


n g 


{x 2 + y 2 + ~z 2 f/ 2 


4A g 


y 


tt g 


{x 2 + y 2 + z 2 f/ 2 


AA g 


z 



Tj %Tj VL 2 g rj' 

Then equations (IT5|) can be simplified to the following dimensionless form 

~ x » _ 2 y> - ^L~x = - 

y" + 2x' = ~ 

v 2 g _ 

Z + W g Z = ~ Q g (i 2 + y 2 +~z 2 fl 2 ' ^ 25 ^ 

Note that the dimensionless equations do not depend on the specific values of the masses 
Mi and M 2 . Thus the result applies to binaries of any masses. The dimensionless form of 
the zero-velocity surface projected to the x-y plane is 

^ 2+ 7^ = "^pk' (26) 

The zero- velocity contours are shown in Figure [TJ Binaries with Ej < E c in region A have 
bounded motion in that they can never escape from A; we call these bound binaries. All 
others are called escaped binaries. 
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-3-2-10123 

x/r . 



Fig. 1. — : Zero- velocity contours in the plane parallel to the Galactic disk, as given by 
equation (|26|) . The red line is the critical contour on which the effective potential $ c g = E c . 
In regions A, B and C, $ e ff < E c and in regions D and E, $ e g > E c . Binary stars with 
large separations can either have $ fr < E c or $ e fj > E c , depending on their positions on 
this plot. We define bound binaries as those in region A with Ej < E c and call all others 
escaped binaries. 
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Once we know the velocity of the center of mass r cm and relative velocity r, we can 
calculate the velocity of each star in the rotating frame from the relations 

r 1 -r cm+Mi + M2 r, r 2 - r cm - _ —v. (27) 



2.2. Kicks from other stars 



In order to study the evolution of the binary systems, we must include the effect of 
encounters with passing stars and other perturbers (e.g., molecular clouds). In this paper 
we only discuss the effects of encounters with stars, but we return briefly to the effects of 
molecular clouds in the discussion of 

As the velocity dispersion of the perturbers (~ 30 km s _1 ) is much larger than the 
velocity difference of the two stars in the binary system (< 1 km s^ 1 ), we can use the impulse 
approximation, i.e., the encounter with the perturber provides an impulsive kick that changes 
only the velocity, not the position, of the subject star. For computational efficiency, we do 
not follow individual encounters but instead consider the total effect of the encounters on the 
binary system after some time interval At p , which is generally large enough to include many 
encounters (see §3.11 for further discussion of this approximation). Let the change of velocity 
of the subject star after this time interval be Av. According to the central limit theorem, 
the effect of a large number of kicks will be the same as that of a Gaussian distribution with 
the same mean fi = (Av) and covariance matrix C a p = (Av a Avj3), where the subscripts 
a, (3 refer to the x, y, z directions. The values of \i a and C Q/ 3 after the time interval At p 
can be computed from the diffusion coefficients 

/i a = D[Av a ]At p , C a(3 = D[Av a Avp]At p . (28) 



We assume that the number density of perturbers with mass in the range M p — > M p + 
dM p is n(M p )dM p , and that the velocity distribution of the perturbers relative to the center 
of mass of the binary is isotropic and Maxwellian, 

dn = f(v p )dM p dv p = ^^ e - v p^dM p dv p , (29) 

where a is the relative velocity dis persion. Expressions for the diffusion coefficients are given 



in equations (7.89) and (7.92) of iBinney fc Tremaind (120081 ). The actual relative velocity 



distribution is more complicated, both because the distribution of stellar velocities in the 
solar neighborhood is triaxial and because the center of mass of the subject binary star has 
its own epicyclic motion, but we do not believe that these complications will alter our results 
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significantly. If the velocity of subject star i relative to the center of mass of the binary is 
Vj, then 



D[Av iia ] =^D[Av {l } t , 

Vi 

D[Av i!a Av itP ] =^£{D[(A Vll )% - l-D[(Av±) 2 ]i} + ^^[(A^) 2 ],, (30) 

where D[Aui|] is the mean change of velocity per unit time along the velocity vector direction 
v, while -D[(At>_i_) 2 ] and D[(Aw||) 2 ] are the mean-square changes per unit time in the velocity 
perpendicular and parallel to v. 



In the limit |vj| <C a, we have 

D[Av\\]i 



D[(Av\ 



\2l 



A\f2TxG 2 (M i p 1 + p 2 ) In A 

: n v * 

3a d 

8 v / 2^G 2 p 2 In A 



D[(Av ± )% = V - P2 , (31) 



where 



Here A is defined to be 



3a 

91 
3a 

Pk = I n{M p )M k p dM p . (32) 



2 

typ 



A = ty l , (33) 

G{M l + M P Y 

where 6 max is the maximum impact parameter considered, v typ is the typical relative velocity 
and M p is the typical perturber mass. Since these parameters enter only logarithmically, we 
can just assume M p = M® and v typ — a for simplicity, and the maximum impact parameter 
frmax can be chosen to be the half of the separation of the two stars when we apply the kic 1 ^ 



By the central limit theorem the distribution function for Av is 

(34) 



/(Av)= ( 27 r)3i|C|V2 ex P 



-i(Av-/i) T -M-(Av-/u) 



Here, Av and \i are 3x1 matrices and M = C _1 is a 3 x 3 matrix. Then the evolution 
of the binary system is followed numerically by repeating the following steps: (i) follow the 
orbital evolution for a time interval At p using the equations of motion (|25|) ; (ii) for each of 
the two stars, draw a random kick velocity Av from the distribution (|34"|) and add this kick 
velocity to the velocity of the star. 



3 We have checked that even though the separations of binary stars have a wide range, different choices 
of 6 max will not change the results significantly. 
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2.3. The diffusion approximation for small semi-major axes 



When the semi-major axis a of the binary system is small enough (a <C rj), the effect 
of the Galactic tide is small compared with the mutual gravitational force of the two stars. 
Then the binary evolves as an isolated two-body system subject to kicks from other stars. 
Moreover the energy kicks from passing stars are small compared to the binding energy of the 
binary, so the evolution can be treated usi ng the diffusion approximation. This problem has 



been studied by previous researchers (e.g.. lKindll977l : iRetterer fc Kinglll982l ; IWeinberg et al. 



19871 ). so we just give the equations here. We will use this diffusion approximation both to 
speed up calculations of the binary evolution at small semi-major axis and to provide insight 
into the numerical results. 

The energy of the binary system E is related to the semi-major axis a by E = —G(M\ + 
M2) / (2a). We define n(E, t)dE to be the number of b inary systems with en ergy in the range 
[E, E + dE] at time t. The diffusion equation reads ([Weinberg et al .1119871 . eq. Bl) 



dn(E,t) 
dt 



dn(E,t) 2 d 2 



dE 



[En{E,t)\ 



where 



1 



SttG z P2 { — )lnA . 

Vrel . 



(35) 



(36) 



Here (1/V re \) is the average inverse relative velocity between the binary and the perturbers, 
which is a/2/7to" _1 under the assumption that the relative velocity distribution is given 
by ( 1291) . As e depends on the binary energy very weakly (through A), we take it to be 
independent of E. We define two dimensionless variables r and h by 



T 



let 
3 1 .Ed 



E 



\Ei\ 



(37) 



where E\ is a scaling parameter. Then h(E, t)dE = h(h, r)dh. With the boundary condition 
n(Ei,t) = an d the initial condition h(E, 0) = 8(E — Eq), the solution to the diffusion 
equation ( l35i) is (IWeinberg et al.lll987l . eq. B13) 



n(h, t) = 512irh 



5/2 



dk 



e- kT k 5 F k (h )F k (h) 
9 + 12k + lQk 2 



(38) 



where h = —E /\Ei\ and 

-5/2 

F k (h) ' ' 



2Vkh 



J5/2 hVkh) J_5/2 [2y/k) - J5/2 (2y/k) J-5/2 hVkh 



(39) 
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with J v {z) the Bessel function of order v. Then the probability that the energy of the binary 
system is larger than E\ for the first time in the interval (t, t + dt) is p(t)dt, where 

p(t) = -— dhn(h,r) = dk ^~i+ . (40) 

W dtj 1 V ' ; 3|Ei| J 9 + 12k + 16k 2 v ; 

To use these results to accelerate our calculation, we choose two semi-major axes ao and ax, 
which are small enough (a < a\ <C rj) that the influence of the Galactic tide is negligible. 
We set E and E\ to be the corresponding binary energies. Typically a = 0.5a! in our 
simulation. If the semi-major axis of the binary system in our Monte Carlo simulation 
random walks to a value smaller than ao at time to, then we draw a random time t from the 
distribution function (jlOjl . which is a fair sample of the time the binary system needs to go 
from ao to a\. If to + t is smaller than the total time of our simulation (10 Gyr), then we just 
give the binary system semi-major axis a\ and other randomly chosen orbital elements as 
described at the start of the following section, and continue to evolve the binary numerically 
from time to + t. If its semi-major axis becomes less than a a second time, we just repeat the 
above calculation. If to + 1 is larger than the total time of the simulation, then we conclude 
that at the end of our simulation, the semi-major axis of the binary system is still smaller 
than a\. Then the probability for the binary system to have dimensionless energy [h, h + dh] 
at time t is 

MKr) dh= r+ T^ d " - (41) 
L n{n,T)dn 

Here the time t (and thus r) is fixed by the condition that to + t is the time at the end of 
the simulation. We draw a random number from this distribution and use this to determine 
the semi-major axis of the binary system at the end of the simulation. We include this 
binary system in the final statistical result after we assign random values to the other orbital 
elements of the orbit as described at the start of the following section. 



3. Numerical simulation of the binary systems 

We simulate the evolution of binary systems for up to 10 Gyr under the influence of 
the Galactic tide and kicks from passing stars. To determine the initial relative position 
and velocity of the two stars in the binary system, we first choose a semi-major axis a, as 
described below. We then choose the inclination angle 6 between the plane of the orbit and 
the Galactic plane randomly so that cos# is uniformly distributed between —1 and 1, which 
corresponds to a spherical distribution. We choose the eccentricity e of the initial orbit so 
that e 2 is distributed uniformly random between and 1, which corresponds to an ergodic 
distribution on the energy surface. The angle between the projected major axis of the orbit 
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on the x-y plane and the x-axis is uniformly distributed between and 2ir. The initial phase 
or mean anomaly of the orbit is also uniformly distributed between and 2n. 

We carry out six simulations of N = 50, 000 binary stars each. In the first four simu- 
lations the systems are "formed" at initial times to that are uniformly distributed between 
and 10 Gyr and followed to t = 10 Gyr, to represent the current state of a population 
with uniform star-formation rate. The initial semi-major axes in units of the Jacobi radius 
are aj/rj = 0.01,0.05,0.10,0.20. In the final two simulations the logarithms of the initial 
semi-major axes are uniformly distributed (Opik's law, eq. [2]) between log 10 (0.001rj) and 
log 10 (0.5rj); in the first of these the binary formation times to are uniformly distributed 
between and 10 Gyr, while in the second the binaries are all formed at to = 0. We label 
these "Opik 1" and "Opik 2" . In these cases a small fraction (< 7%) of the initial binaries 
have already escaped in that Ej > E c . 

Initially, the center of mass is on a circular orbit, which means r cm = r cm = 0. 

We solve the equations of motion (1251) numerically over the time interval At p us- 
i ng an adaptive fourth-o rder Runge-Kutta method and Kustaanheimo-Stiefel regularization 



( IStiefel &: Scheifeld Il97ll ) . The evolution of the center of mass over this interval is given by 
the solution f|T7|) . At the end of this interval, we know the velocities and positions of the 
two stars. Then we generate random velocity kicks Av; for each star from the distribution 
function (1541) . We add Avj to each star while keeping the positions unchanged. With the 
new velocities and positions as initial conditions, we let the binary system evolve for another 
time interval At p . If the semi- major axis becomes smaller than a specified value do, we 
switch to the diffusion approximation as described in §2.31 until either (i) we reach 10 Gyr 
and stop, or (ii) the semi-major axis exceeds a\ > a , at which point we return to a numerical 
simulation If the initial semi-major axis is smaller than a we start with the diffusion 
approximation. We follow each binary system in this way to the time 10 Gyr. If we are 
using the diffusion approximation at 10 Gyr, we draw a random semi-major axis from the 
probability distribution fj4T|) and assign the other orbital elements at random as described 
at the beginning of this section. 



4 For the simulations with ai/rj — 0.05,0.1,0.2, we chose a^/rj — 0.04,0.05,0.1 and a\ = 2ao- In the 
simulation with ai/rj = 0.01, we initially followed the evolution of all stars using the diffusion approximation 
and switched to the Monte Carlo simulation when the semi-major axis exceeded a\/rj — 0.08. In the Opik 
1 and Opik 2 simulations our procedure depended on the initial semi-major axis: for cn/rj > 0.2 we did 
not use the diffusion approximation at all; for 0.08 < at/rj < 0.2, we used a\ = 2ao = a%\ for ai/rj < 0.08, 
we initially followed the evolution using the diffusion approximation and switched to the simulation at 
ai / rj = 0.08. 
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3.1. Values of the parameters in the numerical simulation 



In this subsection, we give the values of the parameters we chose in the simulation. As 
we focus on binaries in the solar neighborhood, the angular speed fl g , vertical frequency 
Vg, Oort constant A g , and epicycle fr equency n g are chosen to b e the values in the solar 
neighborhood, taken from Table 1.2 of iBinney &: Tremaind (120081 ) 



n 



9 

9 



236 km s _1 /8kpc = 9.56 x 10~ 16 s' 1 
2.3 x lO^s" 1 , 



A g =14.8 km s^kpc" 1 = 4.796 x 10 ' s 



37km s _1 kpc _1 = 1.2 x 10 ' s 
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16 „-l 

-1 



(42) 



With these units, the Jacobi radius 



rj = 1.70 pc 



M 1 + M 2 

2M r . 







1/3 



The corresponding velocity and acceleration are 

'Mi + M 2 ^ 1/3 



o TV = 0.050 km s" 1 



2M, 



© 



n 2 g rj = 4.8 x 10 
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km s 



M l + M 2 
2M 



1/3 



(43) 



(44) 



The typical one-dim ensional velocity dispersion in the solar neighborhood is 29 km s _1 
( jDehnen fc Binneyill998l ) and we choose the relative velocity dispersion to be y/2 times this, 
so a = 40 km s _1 . 

The mass function n(M p ) of the stars in the solar neighborhood is given in equation (1) 
of iKroupa et al.l (119931 ): 







n{M p ) = n < 



(M p /M pfi y ai 
(M p /M pfl )~^ 
( {M p , r /M p ^{M p /M Ptr )~^ 



The parameters in this equation are 



if M p <M ph 
if M p j <M P < M Pt0 , 
if M pfl <M P < M P)T} 
if M p>r < M p . 



(45) 



a\ =1.3, «2 = 2.2, «3 = 4.5, 
n =0.087 pc^M" 1 , 



M P)l =0.07M fl 



0; 



M, 



p,0 



O.5M , M. 



p.r 



1M . 



(46) 



-16- 



The moments of the mass function are then 

po = J n(M p )dM p = 0.14 pc~ 3 , 

Pi=J n(M p )M p dMp = O.O45M pc~ 3 , 

p 2 = J n{Mp)M 2 p dMp = 0.029M2 pc ~ 3 . (47) 

Although the dimensionless equations of motion (I25j) and the diffusion coefficients 
_D[(At>||) 2 ] and D[(Av±) 2 } do not depend on the specific values of the binary component 
masses Mi and M 2 , the diffusion coefficient £)[Au||] (eq. [31]) actually depends on these 
masses. When we calculate this kick we choose Mi = M 2 = M & . 

We now describe the choice of the interval At p between kicks. If the two stars are bound, 
we can find the semi-major axis a from the energy equation 

= G(M 1 + M 2 ) = 1 _ G(M 1 + M 2 ) 

2a 2 r 1 ' 

Then the orbital frequency of the binary system Qj, is 



^jm+w, (49) 



a" 



and the orbital period is 



27r _ „r / a 



3/2 

n = ^ = 3xl0 6 yr(— — ) . (50) 
i2 b VO.lpc' 



For a typical number density of stars in the solar neighborhood 0.05 pc -3 and a typical 
relative velocity of the stars 40 km s -1 , the collision time (time between encounters with 
impact parameter less than a) between the binary system and the field star is 

t coll = 1.25 x 10 7 ( ^22) JT . (51) 



a 

If the energy E is positive, the time interval At p is chosen to be 

At p = — = 3.3 x 10 6 yr- (52) 

If the energy is negative and the collision time is longer than the orbital period, the time 
interval At p is just chosen to be the collision time. If the collision time is shorter than the 
period, At p is chosen to be 

0.1 , . 

= ^— prv • (53) 

max(i2 5 , lib) 
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We assumed in §2.21 that the interval At p was large compared to the encounter time. This 
assumption is not correct for bound binaries with semi-major axes < 0.5 pc. Nevertheless, 
our results should accurately reproduce the evolution of the binary so long as the evolution 
time is much longer than the encounter time, since the central limit theorem implies that 
any distribution of velocity kicks with the correct mean and covariance matrix should lead 
to the same cumulative effects. We have checked this by varying the value of At p by a factor 
of 3, and found almost no change in the final distribution of the binary systems. 

We have also checked our simulation code by following 5000 binary systems with initial 
semi-major axis 0.05rj, stopping each simulation when the semi-major axis reaches O.lrj. In 
this range of semi-major axes the Galactic tidal force is at least 10 3 times smaller than the 
gravitational force between the binary components, so the diffusion approximation given in 
§2.31 should be quite accurate. We compared the cumulative distribution of stopping times 
to the distribution predicted by the diffusion approximation (eq. |4Q|) and the maximum 
difference was only 2%. 

3.2. Results from the numerical simulation 

The spatial distributions of the binary stars after 10 Gyr are shown in Figure [3 In each 
panel the relative position r = r x — r 2 is projected onto the x-y and y—z plane. We can see 
tidal tails along the y direction (i.e., the direction of the binary's Galactocentric orbit), which 
can extend to several thousands of Jacobi radii. As the initial semi-major axis a, increases, 
the fraction of stars found in the tidal tails and the maximum extent of the tidal tails both 
grow. In the case a, = O.Olrj, only 11 binaries of the original 50,000 have separation greater 
than lOoj = O.lrj and only 2 have separation greater than rj. In contrast, when a, = 0.2rj, 
70% of the binaries have separation greater than lOa, = 2rj after 10 Gyr. 

The distributions of projected separation (as viewed from a randomly chosen position) 
for the binary stars in the six simulations are shown in the histograms of Figure [3J The 
blue histograms show the full sample while the red histograms show binaries with Jacobi 
constant Ej > E c at the end of the simulation. The figures also show the initial distributions 
of separations in green. Binaries at large separations with Ej < E c must be in regions B or 
C of Figure [TJ while binaries at large separations with Ej > E c may be in any of regions B, 
C, D, or E. 

Remarkably, rather than a cutoff in the distribution of binaries at large separations, we 
see a local minimum in the density, at a projected separation of about 5rj. (For binaries 
with initial semi-major axis do = O.Olrj, the minimum is poorly defined as there are only 
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Fig. 2. — : The relative position r x — r 2 in 50,000 binary systems at time 10 Gyr with 
four different initial semi-major axes: (a) O.Olrj = 0.017 pc, (b) 0.05rj = 0.085 pc, (c) 
O.lrj = 0.17pc, (d) 0.2rj = 0.34 pc. In these simulations the binaries are formed at a 
uniform rate between t = and t — 10 Gyr. We also show two simulations in which the 
initial semi-major axes are uniformly distributed in the log between O.OOlrj = 0.0017 pc 
and 0.5rj = 0.85 pc: (e) formation at a uniform rate between t — and t = 10 Gyr; (f) all 
binaries formed at t — 0. The frames labeled x and y are distributions projected onto the 
x-y plane (parallel to the Galactic plane) at various scales while the frames labeled z and 
y are distributions projected onto the z-y plane. The larger the initial semi-major axis, the 
more stars are found in the tidal tails. 
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Fig. 3. — : Histograms of the projected separation r p of the binary stars after 10 Gyr, for 
initial semi-major axes a« = O.Olrj, 0.05rj, O.lrj, 0.2rj (top four panels) and two cases in 
which log 10 cij is uniformly distributed between log 10 (0.001rj) and log 10 (0.5rj): the initial 
time t is uniformly distributed between and 10 Gyr (Opik 1; bottom left panel), and the 
initial time t is fixed to be (Opik 2; bottom right panel). The projected separation is 
obtained by assuming that the line between the two stars has a random angle to the line 
of sight. The histograms for real three-dimensional separations r are very similar to the 
histograms of projected separations shown here. The blue histograms show the total sample 
of 50, 000 binary stars while the red histograms show the stars with Ej > E c at time 10 Gyr. 
The initial distribution is shown in green. There is a minimum in the distribution near 5rj 
in each case. 
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Fig. 4. — : Distribution of separations of escaped binary stars at 10 Gyr (the precise definition 
of "escape" is given in the caption of Figured]). The figure shows the probability density 
of the escaped stars for the simulations with initial semi-major axis aj = 0.05rj, O.lrj, and 
0.2rj, normalized so that the integral of the probability density over log 10 (r/rj) is unity. 
There are two peaks separated by a minimum around 5rj. The larger the initial semi-major 
axis, the larger is the amplitude and centroid of the exterior peak. The simulation for initial 
semi-major axis aj = O.Olrj is not shown because the number of escaped stars is too small. 
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Fig. 5. — : Escape age At versus separation of the binary system r at the end of the sim- 
ulation. The initial conditions for each panel are the same as in Figure [31 The escape age 
is the interval between the first instant when the Jacobi constant of the binary system is 
larger than E c and the end of the integration at 10 Gyr. The black points are binary stars 
with Ej < E c at the end of the simulation while the red points are binaries with Ej > E c . 
The black points with At > are binary stars that have diffused back to Ej < E c even 
though their Jacobi constants Ej were larger than E c at the time they escaped. The black 
points with At = are stars that never escaped. There is a small concentration of points in 
the lower right panel at Q g At ~ 300; these represent binaries that were in escaped orbits at 
birth. 
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Fig. 6. — : Relative velocity and separation of the binaries at 10 Gyr. The initial conditions 
for each panel are the same as in Figure [31 The red dotted line in each picture is the zero- 
energy line for Keplerian systems, v 2 /2 = G(Mi + M 2 )/r or r/rj = 8(A g /Q g )(Q g rj/v) 2 . 
The binary stars with separation much smaller than tj follow the line very well. The green 
points with error bars show the mean and standard deviation in log velocity for various radius 
bins. The escaped binary stars mostly lie within a small velocity range (\Av | < 10Q g rj ~ 
0.5km s -1 ), which provides us with a good method to find escaped pairs. 
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a few escapers.) The distribution shows two peaks on either side of the minimum; we call 
these the "interior" peak and the "exterior" peak. Most of the stars in the interior peak 
are bound (in the sense that they are found in region A of Figure [T] and have Ej < E c so 
that in the absence of external perturbations they must remain in region A forever). The 
fraction of stars in the exterior peak grows as the initial semi-major axis becomes larger 
or the age of the binaries grows. Note that the minimum is present in plots like Figure [3] 
that show number per unit logarithmic separation; plots of number per unit separation are 
approximately flat between 5rj and a few hundred rj but do not show a minimum. Binary 
stars inside the interior peak in Figure [3] roughly follow the initial distributions, shown in 
green. 

In the Appendix we describe a simple analytic model for the distribution of separations 
that fits the simulations reasonably well. 

In Figure HJ we plot the distribution of separations of the escaped binary stars only — 
stars that are outside region A of Figure [1] or inside region A but with Jacobi constant 
Ej > E c at 10 Gyr — for three different initial semi-major axes a = 0.05rj, O.lrj, 0.2rj. As 
in Figure [3] there is an "exterior" peak, and both the height of this peak and the separation of 
the centroid of the peak grow with the initial semi-major axis of the binaries. More surprising 
is that the distribution of escaped stars in Figure 0] also exhibits an "interior" peak centered 
at r ~ 0.5rj. The orbits of the stars in this peak resemble those of the retrograde irregular 
satellites of the giant planets, most of which are also formally "esca ped" in the sense that 



their Jacobi constant Ej > E c (jHenonlll970l ; IShen &: Tremaind 120081 ) . but nevertheless can 



remain within r j for very long times. Integrations for an additional 40 Gyr, in which kicks 
from passing stars were turned off, showed that the number of stars in the interior peak 
declined with time only slowly, as t~ al . 

We expect that the sooner the binary is disrupted — in the sense that kicks from passing 
stars cause the Jacobi constant to random walk to a value exceeding E c — the larger the 
separation of the binary system will be at the time 10 Gyr. We label the interval since the 
Jacobi constant of the binary first exceeded E c until 10 Gyr (the "escape age") as At (if the 
Jacobi constant never exceeds E c we set At = 0). The relation between the separation at 
10 Gyr and the escape age is shown in Figure The red points have Ej > E c at 10 Gyr 
while the black points have Ej < E c . As noted earlier in Figures [3] and HJ there is a gap 
around the separation 5rj in each panel. The black points with At > had Ej > E c at 
some point in their history, but subsequent perturbations kicked them back to Ej < E c 
(black points with r > Tj must lie in regions B or C in Figured]). The black points along 
the axis At = never escaped, i.e., Ej < E c for the entire integration. For the binary stars 
with separation r ^> rj, the general trend is that the sepration grows with At. The upper 
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Fig. 7. — : RMS line-of-sight relative velocity of the binaries as a function of projected 
separation, at the end of the simulations. The horizontal axis is the projected separation 
normal to a randomly chosen line of sight, while the vertical axis is the RMS line-of-sight 
relative velocity in each separation bin. In Keplerian motion we expect (v^) 1 ^ 2 oc rp 1 ^ 2 , 
shown by the straight line. The relation between the line-of-sight relative velocity and the 
projected separation deviates from the Keplerian relation for r p >rj. 
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Fig. 8. — : Phase-space density of simulated binary stars compared with field stars in the 
solar neighborhood. The binary stars in each simulation at 10 Gyr are divided into different 
bins according to separation. The vertical axis is an indicative phase-space density, defined 
as the spatial density of binaries in each radius bin divided by (t> 2 ) 3 / 2 , the cube of the 
RMS relative velocity. The density is normalized to the total number of binaries (50,000) 
in each simulation, and thus represents the phase-space density that would be observed in 
a catalog of 10 5 / / stars where / is the fraction of stars in the catalog that belong to wide 
binary systems. The blue horizontal line is the analogous indicative phase-space density for 
the field stars, computed as p / (v 2 ) 3 ^ 2 where p = 0.68/Vj is given by equation (J171) and 
(v 2 ) = 3cr 2 with a = 40 km s _1 as derived in §3.11 The phase-space density of binaries 
exceeds the density of field stars out to separations of ~ 10 2 rj. 
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envelope of the points in Figure[5]is roughly r oc (At) Q with a = 1.4-1.5. This behavior has a 
simple physical explanation: the relative velocity random walks due to stellar perturbations 
and therefore grows as v oc (At) 1 / 2 , so the separation grows as r ~ vAt oc (At) 1 ' 5 . 

The relative velocity v and separation r of the binary systems at the end of the simulation 
are shown in Figure [6j The red dotted line is the zero-energy line for Keplerian orbits, 
v 2 /2 = G(M\ + M 2 )/r. In the top four panels, binary stars with separation r much less than 
the initial semi-major axis a; follow this line quite closely, since they are generally found at 
r <C a, only when they are near the pericenter of near-parabolic orbits. In the bottom two 
panels, the binaries follow the zero-energy line closely when r <C O.OOlrj, the lower cutoff 
to the semi-major axis range in the assumed initial Opik distribution. As the separation 
increases, to ~ rj, the typical velocity decreases but the logarithmic spread in velocities 
grows, as shown by the green error bars. 

In Figure [7] we plot the relation between the RMS line-of-sight relative velocity and 
projected separation, as seen from an observer with a random orientation. Different initial 
semi-major axes yield almost the same curve. When the separation is < rj, the relative 

— 1/2 

velocity decreases with increasing separation as r p , as one would expect for Keplerian 
motion. When the separation is larger than rj, the RMS relative velocity increases with 
increasing separation. The minimum RMS line-of-sight relative velocity is ~ Q g rj. 

The maximum relative velocity for separations 3> rj is a few times VL g rj ~ 0.05 km s" 1 
(eq. I44p . two orders of magnitude smaller than the typical relative velocity between unrelated 
stars in the solar neighborhood. Because of this, surveys that provide accurate velocity data 
have far greater ability to identify binaries with r > rj than surveys with only position^. 
To illustrate this, in Figure [8] we plot the indicative phase-space density (number density 
divided by (v 2 ) 3 ^ 2 ) of companions from our simulations, which contain 50,000 binary stars at 
birth. If a fraction / of stars are found in wide binaries, the total number of stars in a catalog 
that is required to obtain 50,000 wide binaries is 1 x 10 5 //. The horizontal line shows the 
analogous indicative phase-space density of field stars in the solar neighborhood, po/(v 2 ) 3 ^ 2 
where po is given by equation (j4"T|) and (v 2 ) = 3a 2 with o = 40 km s -1 as derived in §3.11 
The phase-space density of binaries exceeds the density of field stars out to separations of 
~ 10 2 r j or well over 100 pc. Thus a statistical measurement of the distribution of binaries 
at ~ 100 pc separation can be achieved by a survey such as GAIA that is (i) large enough to 
contain > 10 5 stars that were originally in wide binaries; (ii) accurate enough that the errors 
in distance and velocity are smaller than the separations and relative velocities (~ 100 pc 



5 D. Fabrycky points out that unbound pairs of asteroids, possibly for med by collisional disruption of l arge 
parent asteroids in the past, have been detected by similar techniques (|Vokrouhlickv fc Nesvornvll2008f ). 
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and 0.1-0.2 km s _1 ). 

4. Discussion and conclusions 

We have studied the evolution and disruption of wide binary stars under the gravita- 
tional influence of passing field stars. There have been many treatments of this problem 
already (see the Introduction for references) but most of these (i) ignore the Galactic tidal 
field; (ii) define the binary to be "disrupted" when the Keplerian energy becomes positive or 
when the separation exceeds the Jacobi or tidal radius, and assume that the stars disappear 
instantaneously once they are disrupted. The novel features of our treatment are that we 
include the effects of the Galactic tidal field and follow the evolution of the stars after they 
are disrupted. 

Our simulations show that the usual treatment of binary disruption is oversimplified. 
In particular, 

• The number of binaries does not drop to zero when the separation exceeds the Jacobi 
radius rj\ rather there is a minimum in the density (number per unit log separation) at 
a few times rj, almost independent of the initial semi-major axis distribution. Interior 
to this minimum there is a peak in the density due to the binaries that have not yet 
escaped, and exterior there is a peak due to binaries that are slowly drifting apart 
(Figure ED. 

• Many binaries that have achieved escape energy (more precisely, that have Jacobi 
constants that exceed the critical value E c defined in eq. 1221) remain at separations less 
than the Jacobi radius for many Gyr, either because they are on stable orbits that 
do not escape to infinity or because subsequent perturbations from passing stars bring 
their Jacobi constant back below E c before they have time to escape (Figs. H]and[S]). 

• Because the escaped binary components have small relative velocities, they contribute 
strongly to the phase-space correlation function in the solar neighborhood. Large 
astrometric surveys that can measure three-dimensional distances and velocities to 
sufficient accuracy (~ 100 pc and 0.1-0.2 km s _1 ) can detect this correlation signal out 
to hundreds of parsecs. 



These calculations could be improved in several ways. Our simulations do not include 
perturbations from passing molecular clouds, which a re comparable to the pe r turbations from 
passing stars at a ~ 0.1 pc within the uncertainties (IHut fc Tremaindll985l ; IWeinberg et al. 
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19871 ; iMallada fc Fernandezll200ll ). Moreover the qualitative effects of molecular clouds may 
be different because the impact parameter of the most important cloud encounters is much 
larger than the binary's Jacobi radius, whereas the most important stellar encounters have 
impact parameters smaller than the Jacobi radius. Molecular clouds have a much smaller 
scale height than old stars, so the effects of passing clouds and stars may be disentangled 
observationally by examining variations in th e binary distributi on with the vertical amplitude 
of the center-of-mass motion of the binaries (ISesar et al.l 120081 ) . 



The use of the central limit theorem to model stellar kicks as a Gaussian distribution 
(eq. [M]) is a plausible first approximation but should eventually be replaced by a Monte Carlo 
model of the kicks from individual passing stars. As described in the discussion following 
equation (|53|) the assumption that there are many kicks per interval At p is not correct at 
small semi-major axes. Also, for orbits near the critical Jacobi constant E c there may be 
chaotic phenomena such as re sonance sticking that can o nly be modeled using the actua l 
distribution of velocity kicks (IFukushige fc Heggid I200CH ; iHeggid l200ll ; Ernst et ajj 120081 ). 
Despite these concerns, the tests we have carried out in §3.11 suggest that our results are not 
sensitive to the specific value of At p . 

The distinction between evolution due to the Galactic tidal field ( §2.11) and evolution 
due to impulsive kicks ( § 2 . 2 [) is artificial, since the same stars in the disk contribute both 
the tidal field (apart from a contribution from dark matter) and the kicks (apart from a 
contribution from molecular clouds). The approximation that there is a static tidal field can 
be misleading on timescales less than the collision time (l5"Tj) ; however, we do not b elieve that 
this approximation h as biased our results significantly. See iHeisler fc Tremaind (119861 ) and 
Collins fc Saril (120091 ) for further discussions of this issue. 



There is a large li t erature on tidal tails from star clusters (e.g . . lOdenkirchen et al.ll2001 



Belokurov et al.l 120061 ; iGrillmair fc Dionatod 120061 ; iKiipper et al.l 120081 ) . These differ from 
the binary-star tails discussed here in several ways. Most obviously, clusters contain many 
stars so the structure of the tail from a single cluster can be mapped in great detail; in 
contrast, the tail from a single binary contains only two stars so we must combine many 
binaries to measure the tail properties. A second difference is that the kicks to the orbits 
of stars in a cluster arise from other cluster stars, and therefore cease once the star escapes 
from the cluster, whereas the kicks to a binary arise from passing stars and continue after 
disruption. The most important consequence of this difference is that the length of a cluster 
tidal tail grows oc t, while a binary-star tail grows oc t 1,5 . 

Our results hold only for disk binary stars but it is straightforward to repeat the calcu- 
lation for halo binaries. These are of particular interest because the semi-major axis distri- 
bution of halo binaries can be used to constrain the mass distribution of compact objects in 
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the dark halo flYoo. Chaname fc Gould! 12004 ; buinn et al.1 120091 ) 



We thank Dan Fabrycky, Mario Juric, and Yue Shen for helpful discussions. We also 
thank the referee, Winston Sweatman, for comments that significantly improved the paper. 
This research was supported in part by NASA grant NNX08AH83G, and used computational 
facilities supported by NSF grant AST-0216105. 



A. Diffusion of the escaped binary stars 

Here we give an approximate analytic treatment of our results, by solving for the evo- 
lution of binary stars with r <C r j and r ^> r j separately, then matching the two solutions 
at rj. 

The behavior of binary stars with separation much smaller than r j can be described by 
the diffusion approximation given in §2.31 The probability v e {r)d T for the binary stars to 



escape in the time interval (r, r + rfrjj is (derivative of eq. B19 in [Weinberg et al.lll987l . or 
from eq. HQ] as \E^\ — >• 0) 

»M - V]S^- (A1) 

After the stars have escaped to r ^> rj, the gravitational force between the two stars is 
much smaller than the Galactic tidal force. Then the relative motion in the absence of kicks 
is described by equations (1T91) with the right side set to zero. Moreover their separations are 
dominated by drift along the azimuthal or y direction, as seen from Figure [2], so r ~ y. As 
the amplitude of the epicycle motion is small compared to y, we have y ~ y g , where y g is 
the position of the guiding center of the relative motion (cf. the analogous equations [T7] for 
the motion of the center of mass). Therefore we must determine the equation that governs 
the evolution of y g in the presence of kicks from passing stars. 

The relation between the velocity of the guiding center v g = y g and the relative position 
and velocity (x, y, v x , v y ) is 

v 9 = — A A ' (v y + 2Q g x). (A2) 

This can be verified or derived from the epicycle equations for the relative motion (the 
analogs of eqs. [17] for the center of mass epicycle motion) or from equations (8.101) and 



6 The initial time is set to be zero. 
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(8.102) of iBinney fc Tremaind (120081 ). which relate the orbital parameters and the phase- 
space coordinates to the energy and angular momentum in the epicycle approximation. 

With equation (IA2I and the impulse approximation for the kick, the diffusion coefficient 
for v g is given by 



D[{Av g 



A n 



A n fL 



D[(Av : 



V) 



(A3) 



"9 u "9 - 

The factor of two arises because kicks on both stars contribute to the diffusion of v„ 



Let f(t,y g ,Vg)dy g dvg be the probability that the escaped binary stars lie in the interval 
{Ug) Vg + dy g ) and (v g , v g + dv g ) at time t. Then the distribution function f(t, y g , v g ) satisfies 
the simplified Fokker-Planck equation 



df df 1 r/A , 



1 dv 2 ' 



(A4) 



Here we have neglected the term D[Av g ]df /dv g because df/dv g is small compared to 
d 2 f/dv%. The diffusion coefficient D[(Av y ) 2 } for either star is given by equation (1301) . which 
is nowll 



D\{Av y f] = -iDKAv^} + 



2v 2 



D[{Av ± ) 



(A5) 



For binary stars with large separation, the velocity v is much smaller than a. In this case, 



we have D[(Av±) 2 ] 
(EID and 



2D[(Av\\) 2 ] and thus D[(Av v 



D[(Avu 



D[(Av g 



A n f2 f 
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Then from equations 



(A6) 



'■g J u g , 

Note that the diffusion coefficient is independent of y g and v g so we may label a constant 
D g = D[(Av g ) 2 ]. With the initial condition f(0,y g ,v g ) = 5(y g )8(v g ) and the boundary 
condition that / — > when y g — > oo or v g — > oo, the solution to equation (1A4|) is 



f(t,y 9 , v g 



IXDgt 2 



exp 



^t 3 



+ 



/V 2 



(A7) 



The marginal probability distribution of y g can be gotten by integration over v g , which yields 
f(t>yg) = / dVgf(t,y g ,Vg) = , s I ri ^ ri cx}) 



2-irDgt 3 



3 vi 



2D g t 3 



(A8) 



The subscript i for the diffusion coefficients in equation (|30p is omitted here. 
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Then at the final time £/ = 10 Gyr, the probability that the binary has separation {y g ,y g + 



where p e (t) = ■p e {r)dr j dt is given by ( 1A1I) . 

We compare the probability distribution fl A9j) to the escaped binary stars from our 
simulations in Figure [9] (because the maximum impact parameter 6 max is chosen to be the 
half separation of the binary system at each kick time, which is different at different times, 
we have to choose a "mean" 6 max when we use equation (1A9|) to fit the simulation data). 
We can see that the analytic treatment works quite well at the largest separations, and 
works better if the initial semi- major axis is larger. At small separations the fit is less 
good, presumably because our approximation that the gravitational force between the stars 
is negligible compared to the tidal force is not accurate. 
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Fig. 9. — : Fit of the distribution of the binary stars outside the minimum of the distribution 
in Figure H] (r > 3.16rj) to the analytic model described in the Appendix. All the lines are 
normalized to have unit area. The dashed lines are data from the simulation while the solid 
lines are given by equation (1A9D . The simulation with initial semi-major axis a = O.Olrj is 
not shown here because of the small number of escaped stars. There are two reasons for the 
differences between the simulation and theoretical equation. The first is that we neglect the 
gravitational force within the binary system in equation (jA9ft . which is important near rj. 
The second reason is that the maximum impact parameter 6 max is actually not a constant 
during the simulation while in the theoretical equation we just choose a best fit value of & max , 
which is assumed to be a constant there. We can see the larger the initial semi-major axis 
is, the better our formula can fit the data. 
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